1+1+1 flavor QCD+QED simulation at the physical point 
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We present the results of 1+1+1 flavor QCD+QED simulation at the physical point, in which the 
dynamical quark effects in QED and the up-down quark mass difference are incorporated by the 
reweighting technique. The physical quark masses together with the lattice spacing are determined 
with m n + , m K + , m K o and m n - as physical inputs. Calculations are carried out using a set of 
2+1 flavor QCD configurations near the physical point generated by the non-perturbatively 0(a)- 
improved Wilson quark action and the fwasaki gauge action at /3 = 1.9 on a 32 3 x 64 lattice. We 
evaluate the values of the up, down and strange quark masses individually with non-perturbative 
QCD renormalization. 
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I. INTRODUCTION 

The ultimate goal of lattice QCD is to understand 
the non-perturbative dynamics of strong interactions at 
the physical values of the up, down and strange quark 
masses. Recently we have achieved this goal assuming 
artificial isospin symmetry with degenerate up and down 
quark masses m u = md[l|, I n nature, however, the 
isospin symmetry is broken due to the up-down quark 
mass difference m u ^ and their electric charge differ- 
ence. Their effects are observed in mass splittings among 
isospin multiplets of light hadrons, e.g. m K o — m K ±, 
m n — m p . The magnitude of splittings is tiny yet impor- 
tant since, e.g., it is this difference which guarantees the 
stability of proton. Thus, our next step should be 1+1+1 
QCD+QED lattice simulation at the physical point in- 
corporating the isospin breaking effects. 

A pioneering lattice QCD+QED simulation was car- 
ried out in Refs. employing quenched approximation 
both for QCD and QED, where the photon fields and the 
gluon fields are independently generated by conventional 
Monte Carlo methods and they are superimposed. This 
work was followed by several studies in quenched QED on 
quenched and dynamical SU(3) gauge configurations@- 
flq . These studies have shown that the mass splittings 
among the isospin multiplets of hadrons came out qual- 
itatively consistent with the experimental ones. Never- 
theless, dynamical quark effects in QED should be prop- 
erly evaluated to obtain convincing results. An inclusion 
of dynamical quark effects in QED was first attempted 
in Ref. [ItJ by calculating the ratio of the quark deter- 
minant of QCD+QED to that of QCD on rather small 
lattices using the reweighting technique. Recently results 
for the low-energy constants including dynamical quark 
effects both in QCD and QED are reported [l8|. 

In this article, we present results for 1+1+1 fla- 



vor QCD+QED simulation at the physical point. The 
photon fields are superimposed on the 2+1 flavor 
QCD configurations near the physical point generated 
by the PACS-CS Collaboration|2| employing the non- 
perturbatively 0(a)-improved Wilson quark action and 
the Iwasaki gauge action at j3 = 1.90, corresponding to 
the lattice spacing of a ~ 0.1 fm, on a 32 3 x 64 lat- 
tice. We employ the reweighting technique to incorpo- 
rate the dynamical quark effects in QED[17| and adjust 
the up and down quark masses to their physical values 
independently (2I fig. |20| . We make several improvements 
on the generation of photon fields and evaluation of the 
reweighting factors. 

This paper is organized as follows. Section [TT] is de- 
voted to explain the generation of the photon fields and 
the calculational method of the reweighting factors to in- 
corporate the dynamical quark effects in QED. In Sec. IIIII 
we present the results for 1+1+1 flavor QCD+QED sim- 
ulation at the physical point. Our conclusions are sum- 
marized in Sec. IIVI 



II. REWEIGHTING FACTORS FOR FULL QED 

Given a quenched U(l) gauge configuration superim- 
posed on a 2+1 flavor QCD configuration, the reweight- 
ing factor is evaluated with a set of independent Gaussian 
random noises T)i (i = 1, . . . , N v ) : 
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where -D(e p hQ<j, k*) is the Wilson-Dirac matrix with 
a hopping parameter k* corresponding to the physical 

point and the physical electric charge e p h = ^/47r/137 
multiplied by a fraction Q q specific to each quark fla- 
vor. We employ the determinant breakup technique in 
terms of both the electric charge and the hop p in g pa- 
rameters to reduce the noise fluctuations [U, EH, [HI ■ It 
is clear that the stochastic evaluation of det[W g ] works 
efficiently when the ratio W q = D(e p hQ q , k*)/D(0, K q ) is 
close to the unit matrix, for which local fluctuations of 
the photon field should be suppressed. We employ the 
following procedure to prepare the photon field. 

We first generate the photon field b^(n) on a 64 3 x 
128 lattice, which is twice finer than the QCD lattice of 
32 3 x 64, employing a non-compact pure gauge action in 
the Coulomb g au g e and an appropriate treatment of zero 
modesgSHlil: 

Sphoton = EE -(d,J) v (n) - d u b^(n)) 2 

+EE c '( fl c( 8 AW-w»)) s (3) 



with 

dpKin) = b v (n + fi/2)-b„(n), (4) 

where ft is the unit vector in /z direction on the QCD 
lattice. The coefficient c\ = —0.646 is chosen such that 
the averaged action density on the QCD lattice defined 
by 

i(V M ^(iV)-V^(iV)) 2 \ (5) 
with 

B^(N = n/2) = 6„(n) + &„(n + A/2), (6) 
V^B U {N) = B U (N + A) - B„(N) (7) 

is numerically equal to that of the photon field generated 
with ci = on the QCD lattice. To reduce local fluc- 
tuations we average the photon fields over independent 
paths inside the 2x2x2x2 hypercube on the QED lattice 
embedded in the unit hypercube on the QCD lattice: 
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Finally, we construct the compact link variable J7<j = 
cxp(ie p hQ g -B M ) with e p hQ q the physical electric charge 
of each quark. 

We have checked that the averaging procedure does 
not cause variation of the electric charge by calculating 
the electromagnetic mass splittings among charged and 
neutral pseudoscalar mesons in the flavor non-singlet sec- 
tor: Results obtained with the averaged photon field are 
consistent with those with the conventional ones within 
error bars. We have also found that this procedure siz- 
ably reduces the additive quark mass shift due to the 
QED effects for the Wilson-clover quark action we em- 
ploy. 

Further reduction of the fluctuations in the stochastic 



evaluation of the rewcighting factors is achieved by calcu- 
lating the combined reweighting factor for all three light 
quarks with a single set of noise as expressed in Eqs. ([1]) 
and ([2]) [18]. It should be noted that the odd terms in e for 
l^uds 7 ?*! 2 — 1^1 2 are n °t forbidden by the symmetries in 
the QCD+QED mixed system. The leading contribution 
in the fluctuations is expected to be of order of e. Since 
the total charge for the up, down and strange quarks van- 
ishes, Q u + Qd + Q s = 0, we can expect a cancellation of 
the leading contribution, 0(e) = 0, if the quark masses 
are identical. Even though the degeneracy is broken, the 
partial cancellation of the O(e) terms is expected to re- 
duce the fluctuations in the stochastic evaluation. In or- 
der to illustrate the cancellation we plot the noise depen- 
dence of the odd part in e for IW' 1 ^] 2 — |?/; | 2 (q — u, d, s) 
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FIG. 1: Noise dependence of the odd part in e for \ Wq rfi\ — 
\Vi\ 2 (q ~ u , d, s) with the first breakup interval from e = 
to e p h/400. Total contributions including both the odd and 
even parts in e are also presented (black circles). 



with the first breakup interval e = to e p h/400 in Fig.[TJ 
where the stochastic evaluation is applied to individual 
quarks. The odd part in e is extracted from the calcu- 
lations with ±e. We observe that the fluctuation for the 
up quark is opposite to those for the down and strange 
quarks as a function of the noise, where we expect that 
the next-to-leading contribution of 0(e 3 ) should be neg- 
ligibly small. 



FIG. 2: Configuration dependence of the reweighting factor 
from 2+1 flavor QCD to 1+1+1 flavor QCD+QED at the 
physical point. 



Our results [MeV] Experiment [MeV] 



m K + 



137.7(8.0) 
492.3(4.7) 
497.4(3.7) 
input 



139.57018(35) 
493.677(16) 
497.614(24) 
1672.45(29) 



TABLE I: Hadron masses in physical units. 



III. 1+1+1 FLAVOR QCD+QED SIMULATION 
AT THE PHYSICAL POINT 

A. Parameters 



quark matrix inversions are obtained by a block solver 
algorithm developed in Ref. [23| . which reduce the com- 
putational cost by a factor of 3 to 4 compared with non- 
block solvers. 



Our QCD+QED calculation is based on a set of 
2+1 flavor QCD configurations near the physical point 
generated with the non-perturbativcly 0(a)-improved 
Wilson quark action and the Iwasaki gauge action 
at @ = 1.90 on a 32 3 x 64 lattice by the PACS- 
CS Collaboration @]. The hopping parameters are 
(re ud ,re s ) = (0.13778500,0.13660000). The number of 
configurations is 80 over 2000 MD time units. Using 
the U(l) link variables on a 32 3 x 64 lattice, which 
are constructed as described in Sec. HU 2+1 flavor 
QCD partition function is rcweighted into 1+1+1 fla- 
vor QCD+QED partition function with the hopping pa- 
rameters (re*, re*|, re*) and the physical electric charge e p h- 
The values of the hopping parameters (re*,re*[,K*) = 
(0.13787014,0.13779700,0.13669510) are determined by 
7r + , K + , K° and tt~ masses as the physical inputs. We 
employ 400 determinant breakups for the QED part and 
26 for the fine tuning of the hopping parameters to the 
physical point. Each piece of the divided determinant 
is evaluated stochastically with 12 sets of noises. The 



B. Results 

In Fig. [5] we plot the configuration dependence of the 
reweighting factor from 2+1 flavor QCD to 1+1+1 flavor 
QCD+QED at the physical point, which are normalized 
by the configuration average. The largest value is about 
10 which is the same order as found in Ref. fl8| . 

We calculate the hadron correlators with point and 
smeared sources and a local sink employing the opera- 
tors given in Refs. [l], The smeared source is con- 
structed with an exponential smearing function ^ ( | a?| ) = 
A q exp(-B q \x\) {q = u,d,s) with A uA = 1.2,B u . d = 0.07 
and ^4 S = 1.2, B s = 0.18. To reduce the statistical er- 
rors we employ 16 source points located at every 4 time 
slices with different spatial sites. Figure [3] shows effective 
mass plots for n + ,K + ,K° and fl~ propagators with the 
smeared source. The three horizontal lines are extended 
over the fit ranges, which are [£ m i n , t m ax] = [13, 30] for the 
pscudoscalar mesons and [6, 12] for Q~ baryon, and show 
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TABLE II: Quark masses and their ratios in comparison with our previous work. All results are given with the non-perturbative 
renormalization factors [25|. The first error is statistical and the second comes from the renormalization factor, while the two 
errors were combined in the previous work. 

This work Refs. [2,25] Refs. [1,25] 



in. 



MS 



[MeV] 2.57(26) (07) 

[MeV] 3.68(29) (10) 



mf s [MeV] 83.60(58)(2.23) 86.7(2.3) 87.7(3.1) 

m^f [MeV] 3.12(24)(08) 2.78(27) 3.05(12) 

m u /m d 0.698(51) 

m s /m ud 26.8(2.0) 31.2(2.7) 28.78(40) 



the resulting hadron masses with one standard deviation 
error band. The statistical errors are estimated with the 
jackknife method with a binsize of 100 MD time units. 
The fit results are summarized in Table Q] The scale is de- 
termined as a" 1 = 2.235(11) GeV using the baryon 
mass as input. 

In Fig. |4] we show the hadron mass ratios, m. K + /itiq- , 
™jr+/ m fl- an d rrix" /rriQ- , as a function of the number 
of noises employed for the stochastic evaluation of each 
determinant breakup. The ratios quickly reach a plateau 
around the experimental values (red broken lines) with 
small number of noises. This shows not only that twelve 
noises are sufficient for evaluation of each determinant 
breakup but also that we have achieved a successful tun- 
ing to the physical point. 

Figure [5] shows the ratio of K° to K + propagators 
whose time dependence is expected as 



(jf°(t)jf"(0)) 
(K+(t)K+(0)) 



Z(l- (m K o - m K +)t) , (9) 



where we assume (m^o —m K +)t <C 1. The negative slope 
indicates m^o > rriK+- The fit result with [i m in,imax] = 
[13,30] gives 4.54(1.09) MeV, which is consistent with 
the experimental value 3.937(28) MeV [24| within the 
error bar. 

In Table [TT] we summarize the up, down and strange 
quark masses renormalized at \i = 2 GeV in the con- 
tinuum MS scheme together with the values from our 
previous workfl], 0, [25[. Note that the present results 
employ the non-perturbative renormalization factor ob- 
tained with the Schrodinger functional scheme [25j], and 
for a proper comparison, the previous results for the light 
quark masses in Refs. 1. [|| are updated with the use of 
the non-perturbative renormalization factor in Table 14 
of Rcf. [25|. The first errors are statistical and the sec- 
ond ones are associated with computation of the non- 
perturbative renormalization factors. We neglect the 
QED corrections to the renormalization factor, whose 
contribution would be at a level of 1% or less. A po- 
tentially large finite size effects in QED are discussed in 



Ref. |22j]. The authors in Ref. [12j investigate the finite 
volume effects on the up, down and strange quark masses 
in (1.8 fm) 3 and (2.7 fm) 3 boxes, the latter of which is 
nearly equal to the spatial volume in this work, employ- 
ing chiral perturbation theory including electromagnetic 
interactions. Based on their analyses we estimate the 
systematic errors for our results due to the finite size 
corrections of QED to be -13.50%, +2.48%, -0.07% for 
the up, down, strange quark masses, respectively. 



IV. CONCLUSIONS 

We have presented the results for 1+1+1 flavor 
QCD+QED simulation at the physical point. The dy- 
namical quark effects in QED and the up-down quark 
mass difference are incorporated by the reweighting tech- 
nique. We have made several improvements to re- 
duce the fluctuations in the stochastic evaluation of the 
reweighting factors. The physical values for the up, 
down and strange quark masses are determined by using 
,m K + ,m K o and toq- as input, and quoted in the 
MS scheme at [i = 2GeV. A direct investigation of the 
finite size effects of QED by changing the spatial volume 
is left as a future work. 
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